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Abstract 

A general adaptive approach rooted in stratified sampling (SS) is proposed for sample-based un¬ 
certainty quantification (UQ). To motivate its use in this context the space-filling, orthogonality, 
and projective properties of SS are compared with simple random sampling and Latin hypercube 
sampling (LHS). SS is demonstrated to provide attractive properties for certain classes of prob¬ 
lems. The proposed approach, Refined Stratified Sampling (RSS), capitalizes on these properties 
through an adaptive process that adds samples sequentially by dividing the existing subspaces of 
a stratified design. RSS is proven to reduce variance compared to traditional stratified sample 
extension methods while providing comparable or enhanced variance reduction when compared to 
sample size extension methods for LHS - which do not afford the same degree of flexibility to fa¬ 
cilitate a truly adaptive UQ process. An initial investigation of optimal stratification is presented 
and motivates the potential for major advances in variance reduction through optimally designed 
RSS. Potential paths for extension of the method to high dimension are discussed. Two examples 
are provided. The first involves UQ for a low dimensional function where convergence is evaluated 
analytically. The second presents a study to asses the response variability of a floating structure 
to an underwater shock. 

Keywords: uncertainty quantification, Monte Carlo simulation, stratified sampling, Latin 
hypercube sampling, sample size extension 


1. Introduction 

Monte Carlo methods are used for uncertainty quantification (UQ) in nearly every field of 
engineering and computational science. They are often (rightfully) criticized for their high com¬ 
putational cost - especially for reliability analysis and other applications where extreme response 
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is of interest. Yet, they remain the most robust and effective means of quantifying uncertainty 
in computational analyses. For this reason, there has been much interest in improving the com¬ 
putational efficiency of Monte Carlo methods. As Janssen [lj points out, there are two means of 
accomplishing this computational savings: 1. Improve the convergence rate of the sampling routine; 
and/or 2. Perform a sample set that is optimally small. This work addresses both of these aspects. 

To achieve an improved convergence rate, we focus on a class of variance reduction techniques 
called stratified sampling ei m that includes the popular Latin hypercube sampling (LHS) method 
Hi- Stratified sampling methods operate by subdividing the sample space into smaller regions and 
sampling within these regions. In so doing, the produced samples more effectively fill the sample 
space and therefore reduce the variance of computed statistical estimators. LHS in particular has 
been widely used for uncertainty quantification (e-g. |3 E|) and its properties have been studied 
extensively [3 El El El [13 H]. True stratified sampling, on the other hand, has not received a 
similar level of attention in the UQ literature - largely due to the desirable properties of LHS for 
many applications - although it is widely recognized as an effective means of variance reduction 
E2 E3 E3|. In this work, some properties of stratified sampling methods (including LHS) are 
studied and the case is made for the use of true stratified sampling for certain classes of problems. In 
particular, there are many low-to-moderate dimensional applications where true stratified sampling 
is competitive with or even more effective than LHS. This motivates the need to extend these 
benefits to high dimensional problems and some insights for future work along these lines are 
provided. 

As for minimizing the sample size, there is no generally accepted means to identify an optimally 
small sample set a priori. In this work, we espouse a general adaptive paradigm for UQ (Figure [TJ 
wherein samples are progressively added and analyses conducted until a user-defined convergence 
condition is met. This is not a new concept, but it is one that is perhaps underdeveloped. Early 
work by Gilman m demonstrated this idea for classical Monte Carlo analyses using simple random 
sampling (SRS) and a recent work by Janssen |I] provides a renewed emphasis on its importance. 
Until recently though, most stratified sampling methods required a priori prescription of the sample 
size and precluded sample size extension. Therefore, it was necessary to conservatively oversample 
- otherwise the entire analysis would need to be restarted if the convergence criteria were not met. 
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Figure 1: General adaptive process for UQ. 


Recent developments in LHS enable sample size extension using so-called Hierarchical Latin 
Hypercube Sampling (HLHS) [16] |T7|[18], Replicated Latin Hypercubes (RLHs) [19!, 20f, and related 
methods [21 1 [22lf23] . These methods, while important, do not afford the flexibility necessary to fully 
capitalize on the UQ paradigm promoted in this work. Hierarchical methods, for example, produce 
- at a minimum - a new sample for each existing sample during the extension. Consequently, sample 
size increases exponentially with the number of extensions performed. RLHs grow linearly with 
the number of extensions (n*+i = 2 m) and require each subsequent LHS to be the same size and 
utilize the same stratification as the original. Since there is no refinement of the strata using RLHs, 
there is a tradeoff between loss of desirable sample properties and sample size growth rate. These 
considerations limit the adoption of a truly adaptive UQ process. 

A few adaptive SS methods have been proposed in the past - primarily in the physics community 
focusing on applications in multidimensional integration |2l, [25]. However, these methods have 
slightly different aims and, again, do not possess the desired level of adaptability. To enable 
effective sample size extension from a stratified design, a new methodology - Refined Stratified 
Sampling (RSS) - is proposed that is conceptually simple and adds as many or as few samples as 
desired at each extension. The method operates by dividing existing strata and generating samples 
in the newly created empty strata. The method is proven to unconditionally reduce the variance 
of the sample set when compared with existing methods for sample size extension in stratified 
sampling. R is demonstrated that the method affords convergence rates for statistical estimators 
that are comparable or superior to LHS for certain classes of problems with the added advantage 
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of maximal flexibility in its sample size extension capabilities. Moreover, through an exploration of 
stratum optimality, we suggest that the RSS method provides a convenient and rigorous avenue for 
identifying a stratification of the input space that samples optimally in the output space of interest. 
This is the basis of a parallel effort for applications in reliability analysis [26]. 

The emphasis of this paper is on developing the RSS methodology, which is then demonstrated 
on low-dimensional applications where the benefits of true stratified sampling are observed. Con¬ 
siderations for its extension to high dimensional problems are discussed and motivation provided 
for future research along these lines. Two specific example applications are provided. 

2. Review of Sampling Methods 

Consider a stochastic system described by the relation: 

Y = F(X) (1) 

where the random vector X = {Xi, X 2 ,..., X n } possesses independent components defined over 
the sample space S describing n input random variables with marginal cumulative distribution 
functions (CDFs) DxX')- Correlated random variables are not considered in this work as it is 
common practice to produce a set of uncorrelated random variables from a correlated set using 
methods such as Principal Component Analysis and the Nataf or Rosenblatt transformations. The 
operator F(-) commonly represents a computer simulation such as a finite element model possessing 
strong nonlinearities and/or instabilities such that Y is difficult to assess probabilistically. The 
following presents a brief review of common Monte Carlo methods for randomly sampling X in 
order to perform statistical analysis of Y. 

For notational clarity, subscript i is used to denote a vector component, subscript k is used to 
denote a stratum of the space, and subscript l denotes a specific sample. Additionally, n refers to 
the total number of vector components (dimension), M refers to the number of strata in a design, 
and N refers to the total number of samples. 

2.1. Simple Random Sampling 

Traditional Monte Carlo methods rely on so-called Simple Random Sampling (SRS) or Monte 
Carlo Sampling in which realizations of X, denoted x/; l = 1, • ■ ■ ,N (samples), are generated as 
independent and identically distributed (iid) realizations on S by: 

xu = Dx 1 i (Ui);i = 1,2,...,n (2) 


4 


where £/* are iid uniformly distributed samples on [0,1]. The realizations x are then applied to the 
system y = F(x) and y is statistically evaluated. 

2.2. Stratified Sampling 

Stratified Sampling (SS) divides the sample space S into a collection of M disjoint subsets 
(strata) ftp-] k = 1, 2,..., M with = S and Cl p r)Cl q = 0;p q. Samples x; = [xn, X 12 , ■ ■ ■, xi n \, l 

1.2 , N are generated by randomly drawing (X^fc=i ^k = N) samples within each stratum 
k according to: 

} * = 1,2,..., n (3) 

where Uik are iid uniformly distributed samples on with $ = B Xi (C%) and $ = DxfiCtt) 

and Clk an d C ik denote the lower and upper bounds respectively of the i th vector component of 
stratum For our purposes, stratification is performed directly in the probability space mean¬ 
ing that the strata are defined by prescribing the bounds and {« on the n-dimensional unit 
hyper cube. 

Beyond being disjoint (f2 p n fl q = 0;p q) and Filing the space = S), there are 

no restrictions on the strata definitions. Strata can be defined with different sizes and shapes 
in general and can even be defined a posteriori based on an existing sample set (so-called post 
stratification) [2], The size of the strata in the probability space, p/,, is equal to its probability 
of occurrence pk = With these considerations in mind, we define three general classes of 

stratified designs. 

1. Symmetrically Balanced Stratified Design (SBSD): A stratified design is said to be symmetri¬ 
cally balanced if all strata possess equal probability pi = pj\ Vi, j and the probability space is 
divided equally in all variables of the design. In an SBSD, each stratum is an n-dimensional 
hypercube of equal size. 

2. Asymmetrically Balanced Stratified Design (ABSD): A stratified design is said to be asym¬ 
metrically balanced if all strata possess equal probability p % = pj\ Vi,j but the strata limits 
are assigned arbitrarily. In the simplest case of a ABSD, each stratum is an n-dimensional 
orthotope (or hyperrectangle). However, general ABSDs may possess polyhedral or other 
arbitrarily shaped strata. 
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3. Unbalanced Stratified Design (UBSD): A stratified design is said to be unbalanced if all strata 
do not possess equal probability (i.e. 3 pi,pj : Pi / pj). 

We refer to these terms throughout the paper and samples drawn from these designs are referred to 
as Symmetrically Balanced Stratified Samples (SBSS), Asymmetrically Balanced Stratified Samples 
(ABSS), and Unbalanced Stratified Samples (UBSS) respectively. 

2.3. Latin Hypercube Sampling 

Latin Hypercube Sampling (LHS) divides the range of each vector component X*; i = 1, 2,..., n 
into M disjoint subsets (strata) of equal probability 12 ^; i = 1,2,... ,rr,k = 1,2,..., M. Samples 
of each vector component are drawn from the respective strata according to: 

= D^.(Ui k )-,i = 1,2,... ,n;k = 1,2,..., M (4) 

where Uik are iid uniformly distributed samples on [£/°, £/*] with = — and = The 

samples x/ = [xu, X 12 , ■ ■ ■, xi n \, l = 1, 2,..., N are assembled by (uniformly) randomly grouping the 
terms of the generated vector components. That is, a term xu is generated by randomly selecting 
from the generated components (without replacement) and these terms are grouped to produce 
a sample. This process is repeated M times. 

Because the component samples are randomly paired, an LHS is not unique; there are ( M\) n ~ 1 
possible combinations. Improved LHS algorithms have been developed to determine optimal pair¬ 
ings that either enhance space-filling or reduce spurious correlation (increasing orthogonality). To 
improve space-filling, several Latin hypercube design methods have been developed that mini¬ 
mize the /^-discrepancy ED, maximize the minimum distance between points (‘maximin’ designs) 
|271 ESI ESUaoj, optimize projection properties to ensure samples are evenly spread when projected 
onto a known subspace m, and minimize integrated mean square error and maximize entropy |32j, 
among others. Similarly, numerous efforts have been made to reduce spurious correlations begin¬ 
ning with Iman and Connover [33j and later Florian [33] who rearrange the matrix of samples based 
on a transformation of the rank number matrix. Huntington and Lyrintzis [9j and Vorechovsky and 
Novak |35| utilized iterative optimization methods to reduce spurious correlations while others use 
orthogonal arrays [36j. Further, several authors have developed methods for constructing orthogo¬ 
nal Latin hypercubes that additionally possess enhanced space-filling properties [37, 38]. Many of 
these designs are complex, laborious to implement, and/or computationally intensive. Meanwhile, 
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there is no consensus on whether space-filling or minimal spurious correlations are preferable but 
certainly the two properties are linked and there may be a trade-off in achieving either. 


2-4- Statistical evaluation & variance reduction 

For the purposes of evaluating statistical properties of the different sampling methods, consider 
the general statistical estimator defined by: 

N 

T(yi,...,y N ) = ^2wig(yi) (5) 

1=1 

where yi = F(xi) and x/ denotes a sample generated according to SRS, LHS, or SS, w\ are weights 
attributed to each sample {wi = lr for SRS and LHS but may differ for SS - see below), and g(-) 
is an arbitrary function. If g(y) = y r , then T represents an estimate of the r th moment while 
g(y) = 1 {y < Y}, where 1 {•} denotes the indicator function, specifies estimation of the empirical 
CDF. In keeping with past notational conventions, we denote Tr, Tl, and T$ as the statistical 
estimates produced from SRS, LHS, and SS respectively. 

For SRS, it is well known that the sample variance of the statistical estimator is given by: 

2 

Var[T R ] = ^ (6) 


where a 2 denotes the variance of g(Y). This estimate serves as a benchmark for comparison with 
the considered variance reduction techniques. 

Stratifed sampling does not necessitate an equal probability of occurrence for each sample. 
Consequently, in order to statistically evaluate the response quantity Y generated using stratified 
samples, probabilistic weights are assigned to each sample according to the probability of occurrence 
of the stratum being sampled. That is, for a sample x; E fifc, the sample weight is defined as: 


P[^k] _ Pk_ 
M k ~ M k 


(7) 


where M k is the total number of samples in subject to Y2k=i = N an d Ylt=i P[^k] = L 
Stratified sampling has been proven to unconditionally reduce the variance of statistical estimators 
when compared to SRS. McKay et al.|lj have shown that, for a balanced stratified design (SBSD 
or ABSD): 

1 M 

Var[Ts] = Var[T R ] - — - T ) 2 (8) 

v k =1 
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where /q i: is the mean value of the response evaluated over stratum k, and r is the overall response 
mean. 

Similarly, McKay et al. [4j have shown that the variance of the statistical estimator produced 
from a Latin hypercube design (Tl) is given by : 

Var [Tl] = Var [T R \ + ^ 1 N n^_ l)n “ T X^ “ r ) ( 9 ) 

where t is defined as before, p p is the mean value of LHS cell p defined by the bounds of the 
strata on each of the n marginal distributions for sample l , and R denotes the summation over the 
restricted space of N n (N — l) n pairs (/i p , p q ) of cells having no cell coordinates in common. Notice 
that LHS does not unconditionally reduce the variance of the estimate although it has been shown 
that the variance is reduced for any function F(- ) having finite second moment when N » n [7j. 

n 

Stein |7j has also shown that the closer E(X) is to additive (F(X.) = -Kj(Xj)), the more the 

i=l 

variance of Tl is reduced. 

To the authors’ knowledge, no comprehensive study exists which compares the variance re¬ 
ductions from Eqns. ([8]) and (§. Although such a study is beyond the scope of this paper, it 
is immediately clear that SS and LHS reduce variance through different statistical mechanisms. 
Thus, different classes of problems exist for which each method will afford superior variance reduc¬ 
tion. LHS, for example, has been demonstrated to perform well for applications where F(X.) has 
strong additive components while SS will be shown to perform well for applications with strong 
interactions among the variables of X. 

3. Space-filling, orthogonality, and projective properties of sample designs 

In this section, we study and compare the space-filling, orthogonality, and projection properties 
of the designs considered in the previous section. It is demonstrated that, under certain conditions, 
SS has merits that, to date, have not been exploited and can lead to meaningful improvements in 
sample design. 

3.1. Space-Filling 

Numerous metrics have been developed to quantify the space-filling of a sample design including 
various discrepancy measures and the maximin and minimax distances. Each metric provides 
a slightly different interpretation of space-filling and can therefore yield apparently inconsistent 



results (one design may be “good” by one measure and “poor” by another). It is not our intention 
to discuss the relative merits of these measures in general but the interested reader is referred to the 
works of Fang et al. [39| and Dalbey and Karystinos pfOj for further discussion. Instead, we aim to 
compare SRS, LHS, and SS using a couple of metrics to gain some insight into how these methods 
fill the space. For this purpose, we will utilize the commonly employed wrap-around /^-discrepancy 
[4TJ and a new method derived from the Voronoi decomposition [42] . 

The Voronoi decomposition is defined such that each point in space is associated with the 
nearest sample point. That is, each Voronoi cell is defined as the region of the space V, denoted 
Rk and associated with point Pk, containing all points in V whose distance to Pk is less than or 
equal to the distance, d, to any other point Pj,j / k. That is: 

R k = {x e X|d(x, P k ) < d(x, P fc )Vj / k} (10) 


where, for the purposes of this work d(x, y) is the Euclidean distance. Dividing the space in 
this way facilitates evaluation of different sample designs because the cell sizes provide a direct 
measure of space-filling. In particular, large Voronoi cells correspond to regions of the space that 
are sparsely sampled while small cells correspond to regions of the space that are densely sampled. 
By this measure, an optimal space-filling design produces cells whose sizes have minimal ensemble 
variance. In low dimension, the Voronoi decomposition can be established explicitly. In high 
dimension, however, the Voronoi decomposition must be estimated implicitly. To do so, the space 
is populated with a large number of points (N p ~ 10 5 — 10 s ) and the samples are each attributed to 
their nearest point in the sample set. The size (volume) of each Voronoi cell is therefore determined 
by the proportion of points attributed to each sample, N t : 


V, = ^. 

N n 


( 11 ) 


With this general procedure for estimating cell sizes, we propose that the coefficient of variation 

of the cell sizes (——), referred to as the “V-metric” herein, provides an effective metric of space- 

Pv 

filling that can be applied regardless of sample size or dimension. It is straightforward to interpret 
(low value means that the points occupy approximately the same amount of space) and is visually 
clear when observed in 2D. Consider 100 samples drawn using SRS, LHS, and SBSS. A plot of 
the Voronoi decomposition (Figure [2]) shows that the SBSS produces cells that are more uniform 
in size than LHS and SRS. The inset statistics verify that, indeed, SBSS produces cells whose 
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Figure 2: Voronoi tessellation of the 2D probability space from 100 samples drawn using SRS (left), LHS (middle), 
and SBSS (right). 




Number of Samples 

(b) 


Figure 3: (a) Voronoi cell size coefficient of variation (V-metric) and (b) maximum/minimum cell size for different 
sample designs as a function of sample size for 2D samples. 


size variance is significantly smaller than both LHS and SRS. These relations are valid regardless 
of sample size as demonstrated in Figure [3j which shows the V-metric along with maximum and 
minimum observed cell sizes (from all sample sets) relative to sample size as evaluated by Monte 
Carlo simulation averaging over 10 independent trials. Note also that, by this metric, the space¬ 
filling of LHS degrades while the SBSS and ABSS maintain their space-filling properties as the 
sample size grows. 

To study the dependence on dimension of the space-filling properties for these sample designs, 
we compute the average V-metric for each design from 1,024 samples (repeated 10 times) of dif¬ 
ferent dimension. Figure [4(a) shows that, by this metric, SS provides a more even distribution of 
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Figure 4: Space-filling properties of LHS, LHS with correlation correction (LHS-Corr), SRS, and SS as measured 
using (a) the Voronoi cell metric (‘V-metric’) and (b) the wrap-around L 2 discrepancy. Note that for n = 2, 5,10 SS 
designs are SBSD while for n = 8, 20 SS designs are ABSD. 


samples than either LHS or SRS for low-dimensional random vectors while all sampling methods 
lose effectiveness for higher dimension. Notice that, for LHS, the use of the procedure by Iman and 
Conover [33] (denoted LHS-corr) to reduce spurious correlation has only a small effect on space¬ 
filling although, as we will see in the following section, it significantly improves the orthogonality 
of the sample set. 

The second metric we consider is the wrap-around ^-discrepancy (Dl 2 ) [S] defined through 
the standard discrepancy measure: 


D(X) = 


xnc- 


M 


— Vol(c M ) 


N 


( 12 ) 


with c M C [0, 1] M such that || • || denotes an L 2 norm and the subset c M includes all hyperrectangles 
that can exist in a periodic domain [0, 1] M . Figure 4(b), which plots Dl 2 for different sample designs 
and dimensions, is in contrast with the Voronoi metric in that LHS is shown to provide better space¬ 
filling properties. This is due, as highlighted by Dalbey and Karystinos [ID], to the fact that Dl 2 
is sensitive to differences in low-dimensional subspaces. Since LHS discretizes the ID subspaces 
very evenly, Dl 2 indicates a high degree of space-filling. Meanwhile the V-metric quantifies the 
uniformity over the whole M-dinrensional space. 

The comparison above does not conclude therefore that one sampling method is “better” than 
the other at space-filling. Instead, it serves to highlight that LHS and SS fill the space in different 
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ways. Owing to the hierarchical ordering principle [43] . which states that main effects and low- 
order interactions are usually more important than higher-order effects, the first-order space-filling 
of LHS is desirable for many problems. However, when interaction effects play an important role, 
the higher-order space-filling of SS may be desirable. 


3.2. Orthogonality 

An orthogonal sample is one whose sample correlation matrix is diagonal. That is, the random 
samples are perfectly uncorrelated. In practice, an orthogonal (or near orthogonal) design can be 
difficult to achieve. Latin hypercube designs often produce spurious correlations unless corrected 
or iterated in some way (e.g. [331 M, 0 [35]). Stratified samples however, intuitively possess little 
artificial correlation to begin with. This is demonstrated in the following where we consider two 
measures of orthogonality similar to those considered by Cioppa and Lucas [38]. 

First, we evaluate the spurious correlations among the sample variates for a two dimensional 
sample produced using the various designs and consider their variations from the intended zero 


correlation using Monte Carlo simulation and averaging over 10 independent trials. Figure 5(a) and 


5(a) show the standard deviation and maximum spurious correlation respectively for five different 


sample designs as a function of sample size. Notice that the stratified samples (both symmetric 





(a) (b) (c) 

Figure 5: Evaluation of spurious correlation produced by various 2D sample designs as a function of sample size: (a) 
standard deviation of spurious correlation, (b) maximum spurious correlation, and (c) condition number. 

[SBSS] and asymmetric with 2:1 aspect ratio [ABSS]) produce smaller spurious correlations than 
even the LHS with explicit correlation control (LHS-Corr) while LHS without correlation control 
produces relatively large spurious correlations. Moreover, the relative improvement of a stratified 
design over LHS and LHS-Corr increases with sample size. 

Next, we evaluate the condition number of X 7 X, where X is the N x n design matrix with 
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values in the probability space scaled to the range [—1,1], defined by: 


cond (X T X) = — 


(13) 


where <j>\ and r/> n are the 1 st and n th eigenvalues of X 7 X respectively. The condition number is a 
more general metric of orthogonality than any single correlation value with a value of 1 indicating 
a perfectly orthogonal sample. It is therefore useful for higher dimensional samples. We consider 


first the condition number for the 2D samples above as a function of sample size (Figure 5(c)) 
and observe the same general trend. However, as the dimension grows, stratified designs are less 
effective at achieving orthogonality as demonstrated in Figure [6] which shows the mean condition 
number from 10 independent trials for samples of size N = 256 (Figure [6(a )J and N = 1024 (Figure 
|6(b)|) as a function of dimension for the different sampling methods. As the dimension grows, LHS 




Figure 6: Condition number of the design matrix as a function of dimension for the different sampling methods from 
(a) N = 256 samples and (b) N = 1024 samples. 

methods with active correlation reduction significantly improves the orthogonality of the sample 
while SS produces sample sets with reduced spurious correlations up to moderate-dimensional 
samples without the need for any sample adjustment. 

3.3. Projective properties and variable interactions 

Often, in an uncertainty analysis, certain variables have a strong effect on the response while 
others are relatively insignificant. Global, or variance-based, sensitivity analysis is one means of 
assessing the significance of each variable (and interactions) [33]. But, this significance can usually 
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only be established a posteriori (i.e. after concluding the Monte Carlo study). Consequently, it 
necessary that a sampling method project all variables through the transformation effectively in 
order to capture the effects of the significant variables (even when they are not known). This 
property of a sample design is referred to as the projective property and it is perhaps the greatest 
strength of LHS. Because LHS discretizes each variable finely (highly resolving all of the marginal 
distributions), it projects each variable through the transformation well. 

The projective properties of LHS follow directly from the properties illuminated by Stein [7], 
who showed that LHS has the effect of filtering out the additive components (or main effects) of 
the transformation. More specifically, Stein showed that if the transformation /i(x) is decomposed 
in terms of its main effects fi a (x) and interaction effects r(x) as: 

h(x) = /i 0 (x) +r(x) (14) 


then the variance of the LHS estimator becomes: 


Var(T L ) = r(x) 2 dF(x) + O 


(15) 


The variance associated with the main effects h a (x) is very small, while the variance 

associated with the interactions r(x) is equal to that of a standard Monte Carlo estimate (i.e. 
there is no variance reduction on the interactions). By contrast, SS reduces the variance on the 
main effects and the interactions in equal measure. This leads to the following conclusion. LHS 
has excellent projective properties for individual variables. SS, on the other hand, has moderate 
projective properties for individual variables but much better projective properties for variable 
interactions. This benefit, however, diminishes as the dimension grows as demonstrated by the 
following example. 

Consider two simple transformations. The first is an additive function defined by: 


Y i = -Yl Xi 

n ' 


(16) 


i =1 


where X{ ~ C/(0,1). The second is a multiplicative function with strong variable interactions given 
by: 

n 

Y 2 = Hx i . (17) 


2=1 


where X* ~ U{\ — \/(3), 1 + ^/(3)). Each transformation has been designed to produce E\Y\] = 


E[Y 2 ] = 1. 
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Figure [T] shows the average standard deviation from 1,000 Monte Carlo estimates of E\Y\ ] and 
E7[Y2] from 1024 samples using SRS, LHS, and SS for varying problem dimensions. Note that 
since E[Xi\ ^ 0, the transformation in Eq. © possesses both main effects and interaction effects. 
This can be seen by expressing Xi in terms of zero mean variables Xj as Xj = Xj — 1 such that 
I 2 = ( X\ — 1)(X2 — 1)... (X n — 1). It is these main effects that enable a variance reduction from 


LHS on Eq. (17). As expected, LHS performs exceptionally well for transformation Y\ regardless 



Figure 7: Response standard deviation for (a) an additive transformation and (b) a multiplicative transformation as 
a function of dimension for SRS, LHS, and SS. 


of dimension. However, SS reduces variance considerably over LHS for Y 2 . This variance reduction 
diminishes with dimension but remains more effective than LHS up to n = log 2 (IV) (here n = 10 for 
N = 1024), after which neither method is capable of producing a meaningful variance reduction. 

4. Sample size extension 

The adaptive UQ methodology used in this paper requires the ability to easily add samples to an 
existing set. Sample size extension for SRS is straightforward because samples are iid realizations. 
However, in SS and LHS, sample size extension is not necessarily trivial. For SS, samples are 
traditionally added to existing strata although, as we will show, this is not optimal. Extension of 
LHS meanwhile requires to maintain equal weight for all samples for which a few methodologies 
have been developed recently. 

The first methodologies developed to extend Latin hypercube samples are referred to as Repli¬ 
cated Latin Hypercubes (RLH) |19)[20j. RLH entails performing an additional (independent) LHS 
with identical stratification upon completion of the prior LHS. Using RLH, the minimum achievable 
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sample size extension is limited by the original LHS size (it requires adding N additional samples 
at each refinement) and the sample does not benefit from any additional space refinement. 

Recent works by Tong, [16], Sallaberry et al. HZ!, and Vorechovsky et al. [THj have proposed 
methods commonly referred to as Hierarchical Latin Hypercube Sampling (HLHS). Although these 
methods differ in their details (e.g. the methods in mm enable extension for samples of correlated 
variables), the basic premise is the same. Given an LHS of size N, the strata of each sample 
component are further divided into t + 1 strata (one containing the original sample) - where t is 
referred to as the refinement factor - and the new components are randomly paired as in a typical 
LHS implementation. HLHS introduces N new = tN new samples to the set and in general, through 
r sample size extensions, the size of the sample set grows exponentially as: 

Ntot = N(t + l) r . (18) 

Its space-filling properties and statistical convergence are improved over RLH though because each 
sample size extension involves a division of the sample strata. 

Lastly, an alternate means of sample size “extension” is through a method called Nested Latin 
Hypercube Sampling (NLHS) developed in |21j and [23j. Extension is presented in quotations 
because it is not truly an extension technique. Rather, it works in the opposite manner as HLHS by 
producing one very large LHS sample (larger than would presumably be necessary for the problem 
at hand) and dividing this sample into smaller “nested” Latin hypercubes. It can therefore be used 
as an extension method by considering first a single nest and then considering subsequent nests as 
the sample extensions. 

These developments for adding samples to a Latin hypercube design represent a significant 
development toward achieving the desired adaptive Monte Carlo framework. Each of the methods, 
with the exception of RLH, strictly maintain the character (and all associated properties) of a Latin 
hypercube design as discussed in Section 3. Given the broad appeal of LHS and its many desirable 
features, these methods are attractive for many practical applications. 

5. Refined Stratified Sampling 

The primary drawback of the LHS extension methods is the rapid sample size growth associated 
with sample size extension. The methodology developed in this section, referred to as Refined 
Stratified Sampling (RSS), makes use of the unequal sample weighting allowed by stratified sampling 
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to extend a stratified sample set by a single sample. Furthermore, after each extension, the sample 
maintains the properties of a stratified sample with weights that can be easily computed if the 
sample components are independent and uncorrelated. 

Consider the input random vector X with n independent and uncorrelated components defined 
on the sample space S as described in Section 2. Given an initial stratified sample set of size N 
distributed over M = N strata, the Refined Stratified Sampling methodology proceeds as follows: 

1. Select a stratum to divide according to the following criteria: 

• If a single stratum exists such that Wf. > Wj Vj / k, divide this stratum. 

• If N s strata fi*,; k = 1, 2,. .., N s exist with Wk = rriax(rcj), randomly select the stratum 

i 

Clk to divide with the probability of dividing fife equal to P[D(Cl^)] = —. 

2. Divide stratum fi*. in half according to the following criteria: 

• Compute the unstratified lengths of fl/ c , defined as: 

Aifc = ^-e«;* = l,2,...,n (19) 

where and £■£ are defined as in Section 

• Determine the maximum unstratified length of as A& = max(Ajfc) 

i 

• Divide stratum along the component i* corresponding to A*,. If N c components exists 
such that individual unstratified lengths A,*/. = A j.;i* = 1,2, ...,N C then randomly 
select the component to divide with the probability of dividing component i* equal to 

pmn] = ^ 

3. Keeping all existing samples x^: l = 1,2 ,,N, randomly sample in the newly defined empty 
stratum. Compute new weights for each of the existing samples and the new sample according 
to Eq. 0. 

4. Repeat steps 1 - 3 for every new extension. 

The RSS process is described graphically in the flowchart provided in Figure [8] for a two- 
component random vector starting with N = M = 1 and showing the first four sample-size exten¬ 
sions using RSS. 


2.2 
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Figure 8: Refined Stratified Sampling: Flowchart of sample size extension procedure for two random variables. 
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5.1. Why refine strata? 

The methodology presented herein utilizes a refinement of the strata definitions rather than 
adding samples to existing strata. The benefits of this strategy can be significant - but can also 
be detrimental if not implemented appropriately. The benefits of the method rely on the following 
theorem. 


Theorem 5.1. Let fix denote a stratum of space S and k = 1,..., N ss denote N ss disjoint 
substrata of Cti such that U k=1 u} k = fix and P[oj k } = P[^j\ Vfc, j € 1 ,...,N SS . The variance of 
a statistical estimator computed using N ss samples drawn singularly from u k ', k = 1,... ,N SS will 
always be less than the variance of the same estimator computed using N ss samples drawn from fii. 


Proof Consider a statistical estimator from a stratified design, Tg, as given in Section 2.4 Next, 


consider two different stratified designs wherein all strata are identical except in one region, denoted 
fix- In design one, the region fi i has only one stratum possessing N ss samples. In design two, flj 
is divided into N ss balanced strata, each possessing a single sample. The variances of the statistical 
estimators from these two designs are given by: 


M 


3= 2 3 


2 m 
^[Tsi\ = + 


Ns S M 2 

Var [T s2 ] = J2 Pk a k + 


k =1 


3=2 3 


N, 3 


(20a) 

(20b) 


where Y2k=i Pk = Pi an d the M term summation refers to all points in the stratified design outside 
of region fix. The difference in variance between these estimators is given by: 


Ns, 


2 

Var [T a x] - Var [T s2 ] = -^pW k 


( 21 ) 


k =1 


Pi 

Under the condition that p k = —— (i.e. balanced stratum refinement of fix), the difference reduces 

JNss 


to: 


Var [T si ] - Var [T s2 \ = 



( 22 ) 


It is straightforward to show that the term on the right hand side of Eq. (22) is strictly positive. 
Therefore, variance is reduced using a balanced restratification of fix- □ 
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Figure 9: Identification of the optimal stratum refinement location for a normal distribution. 


5.1.1. Optimal strata refinement for various output distributions 


Theorem 5.1 states that a balanced stratum refinement of any given stratum will always reduce 
the variance of statistical estimates when compared with simply adding samples to the existing 
stratum. Note, however, that Theorem |5.1| does not say that any stratum refinement in general 
will result in reduced variance. In fact, improper stratum refinement can increase the variance. 


Furthermore, Theorem 5.1 does not imply that a balanced stratum refinement produces the optimal 
variance reduction. It often will not. 

Consider a statistical estimate of the expected value of Y. Ts = E [Y], from a stratified design. 
Here, we establish the optimal stratum division (on the basis of variance reduction) for different 
distributions of Y produced from the transformation Y = -F(X). For demonstration purposes, the 
estimate is computed initially from two samples, one in each stratum Hi and Q 2 with p\ = P [fii] = 
0.5 and P 2 = P [^ 2 ] = 0.5 as depicted in Figure [9j We are interested in identifying the optimal 
location to divide the strata such that a single sample can be added while minimizing Var [Ts] ■ In 
each example, stratum SI 2 is divided into two substrata denoted U 21 and UJ 22 possessing probability 
weights poi and P 22 subject to: 


P 21 = P [w 2 i] = zp 2 (23a) 

P 22 = P [^ 22 ] = (1 - z)p 2 (23b) 

where z E [0,1], referred to as the ‘Imbalance Factor,’ is a measure of imbalance of the resulting 
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Table 1: Optimal strata division for uniform, normal, and lognormal output distributions. 


Dist. 

P21 

P22 

°21 

^22 

z 

t/(o,i) 

0.25 

0.25 

5.2e — 3 

5.2e — 3 

0.5 

N( 0,1) 

0.3162 

0.1838 

6.53e - 2 

0.21 

0.632 

L!V(-1.49,1.27) 

0.4426 

0.0574 

0.1165 

6.97 

0.885 


Table 2: Variance of the estimated mean value for different strata refinements and different output distributions. 


Var[T s ] 



Three Samples 

Three Samples 

Three Samples 

Dist. 

Two Samples 

No Refinement 

Bal. Refinement 

Opt. Refinement 

u( 0 , 1 ) 

0.0104 

7.8125e — 3 

5.859e - 3 

5.859e - 3 

N( 0,1) 

0.18175 

0.1363 

0.1083 

0.1045 

L!V(—1.49,1.27) 

0.4138 

0.2073 

0.1696 

0.04667 


stratum division (Figure [9]). Note that z = 0.5 corresponds to balanced stratum refinement. The 
objective of the optimization can be stated as: 

minimize: Var [T s \ = p\o\ + p\\o\\ + P 22 a 22 

Z 

subject to: P 21 = zp 2 (24) 


Pi + P21 + P22 = 1 

We consider three different output distributions: Y ~ Normal(0,1); Y ~ Uniform(0,1), Y 


LogNormal(—1.49,1.27). Optimization as defined by Eq. (24) yields the stratum refinement out¬ 
lined in Table [l] for each distribution. Table [2] meanwhile compares the resulting variances of T$ 
for each distribution considering different strata definitions: 1. Two samples - one from fij and 
one from CI 2 ', 2. Three samples with no refinement - a single sample is added to stratum O 2 ; 3. 
Three samples with balanced refinement - stratum is divided such that P 21 = P 22 = 0.25; 4. 

Three samples with optimal refinement of stratum Q 2 as defined in Table [l] Notice that optimal 
stratification can have a profound influence on variance reduction. 

Figure [lO] shows the variance of the mean estimates for the normal, uniform, and lognormal 
distributions as a function of the imbalance factor z with the optimal value demarcated by an ‘x’. 
The dashed line shows the variance resulting from samples produced without stratum refinement. 
Notice that the balanced refinement is not necessarily the optimal but always reduces variance 
while many refinement strategies (e.g. z < 0.4 for lognormal output) actually increase the variance 
relative to estimates produced without refinement. Analysis of these results shows that the closer 
the stratum Q 2 is to possessing symmetric conditional distribution Dy-{y\y E ^ 2 ), the closer the 
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(a) (b) (c) 


Figure 10: Variance of the estimated mean value from three samples with (a) normal, (b) uniform, and (c) lognormal 
output distributions as a function of stratum imbalance factor. 


balanced refinement strategy becomes to optimal. Conversely, when the conditional distribution 
Dy(y\y e n 2 ) possesses strong asymmetry, the optimal refinement will be imbalanced. In such 
cases, the balanced refinement procedure proposed herein produces diminished return in terms of 
variance reduction. 

The analysis performed in this section assumes stratification in the domain of the output Y. 
In general, the corresponding strata in the input domain of X cannot be identified unless F(X.) is 
known and invertible. Additionally, the distribution of Y is not known in general. Consequently, it 
can be difficult and sometimes impossible to identify the stratum division for X that corresponds 
to the optimal division on Y - particularly if F(X.) is strongly nonlinear, non-monotonic, discontin¬ 
uous, or otherwise complex. Balanced stratum refinement on the input space, on the other hand, 
necessarily corresponds to a balanced stratum refinement on the output space and will always pro¬ 
vide a reduction in variance. For this reason, its use is recommended when stratum optimization 
is not possible. 

Identifying optimal stratum refinement, however, represents a potentially important direction 
for future research. While some similar concepts have been developed for optimally sampling the 
input space with probabilistically weighted draws to match the input probability model (e.g. m) 
the potential capability to use RSS within the adaptive framework promoted herein represents 
an opportunity to define a partitioning of the input space that optimizes the partitioning of the 
output space. This development is one that is potentially very powerful for efficient uncertainty 
analysis. An initial exploration along these lines is undertaken in a parallel effort for applications 
in reliability analysis [26] where strata in the input variable space are divided such that samples 
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are concentrated in the output space in the region of the limit state function. 

5.2. RSS yields unbalanced stratified designs 

Globally, RSS produces stratified designs that are balanced only when N = M = 2 P , are sym¬ 
metrically balanced only when N = M = 2 np for integer values of p, and are otherwise unbalanced 
(assuming an initial sample size of one). In general, producing unbalanced designs is recommended 
only when the output distribution calls for it (see the previous section). However, the output 
distribution is unknown in general and therefore balanced stratification should be used whenever 
optimization is not possible. For this reason, it is recommended that the initial stratified design 
be a balanced one (asymmetric or symmetric). The proposed methodology will then produce un¬ 
balanced designs during most iterations but this is justified by Theorem 1. That said, the analysis 
conducted in the previous section and the variance reduction expression in Eq. ([8]) together imply 
that optimal refinement will not rely on randomly selecting the stratum to divide but rather se¬ 
lecting it based on some established criteria that may optimize variance reduction. One condition 
that could be considered, for example, is to enforce some symmetry conditions in stratum selection 
such that the stratum design does not become too heavily imbalanced. 

5.3. Expected performance and extension to high dimensional applications 

As demonstrated in Section 5.1, the efficiency of the method is dependent on the resulting 
distribution (and more specifically the conditional distributions of the response from within each 
stratum). Given these considerations, the method is expected to converge more rapidly for condi¬ 
tional output distributions that are closer to symmetric. For problems whose output distribution 
is strongly asymmetric and characterized by a dominant peak and heavy tails (i.e. those with high 
skewness and kurtosis), the convergence rate of the RSS method will be reduced due to the sub¬ 
optimality of the stratum refinement. Optimal stratum refinement methods need to be developed 
to fully capitalize on the variance reduction potential of RSS. 

In its present form, the proposed methodology is not well-suited to problems with very high 
dimensional input vectors. As the analysis in Section 3 implies, stratified sampling loses considerable 
effectiveness for dimensions above n = log 2 (IV) and its benefits are most pronounced only for low 
dimensional problems. To extend the RSS method to very high dimensional spaces, it is currently 
being considered to conduct RSS on low-dimensional subspaces with samples generated on these 
subspaces randomly paired as in LHS. Moreover, a major deficiency of the proposed method from a 
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Table 3: Input distributions and resulting analytical moments of the output. [LN = LogNormal(/x, a), U = 
Uniform(a, b), N = Normal(/r, <r)]. 


Dist. 

Xi 


a 

E[Y] 

Var[Y] 

Skew[Y] 

Kurt[Y] 

A 

LN(0,0.01) 

U(0,20) 

N(l,0.1) 

-113.33 

12012.0 

-0.77 

-0.55 

B 

LN(0,0.1) 

U(0,10) 

N(l,0.1) 

-23.37 

621.18 

-0.89 

-0.24 

C 

LN(0,0.1) 

U(0,7) 

N(l,0.1) 

-9.32 

121.97 

-0.97 

-0.09 

D 

LN(0,0.1) 

U(0,6) 

N(l,0.1) 

-5.98 

58.46 

-1.01 

0.002 

E 

LN(0,0.1) 

U(0,5) 

N(l,0.1) 

-3.31 

23.65 

-1.08 

0.17 

F 

LN(0,0.3) 

U(0,5) 

N(l,0.1) 

-3.10 

25.03 

-1.17 

0.81 

G 

LN(0,0.4) 

U(0,5) 

N(l,0.1) 

-2.87 

26.80 

-0.99 

2.08 

H 

LN(0,0.45) 

U(0,5) 

N(l,0.1) 

-2.70 

28.85 

-0.48 

9.42 

I 

LN(0,0.475) 

U(0,5) 

N(l,0.1) 

-2.60 

30.56 

0.09 

23.84 

J 

LN(0,0.5) 

U(0,5) 

N(l,0.1) 

-2.48 

33.05 

1.08 

60.62 


dimensionality perspective is that rectilinear stratification requires stratum division in only a single 
dimension. A generalized RSS algorithm would refine strata of arbitrary shape (i.e. Voronoi cells) 
to allow refinement in multiple dimensions simultaneously - thereby greatly reducing dimensional 
dependence. Such algorithms are the topic of a forthcoming work. 


6. Applications 

6.1. Statistical analysis of a stochastic function 

Consider the cubic polynomial function given by: 

Y = F(X) = XfX 2 - aX t X\ + X x X 2 (25) 


The function possesses three random variables: X±, X 2 . and a. Several sets of input distributions 
are considered to elicit significantly different response distributions as described through analytically 
computed moments in Table [3j The cases cover a vast range of distributions from those approaching 
uniform to those with sharp peaks and heavy tails. Convergence is defined based on the following 
relative error: 


< 7y - Var[Y] 


Var[Y] 


Y Cth 


(26) 


where Var[Y] is the analytically determined variance, <Jy is the variance computed from the gener¬ 
ated samples, and eth = 0.01 (in general, analytical moments will not be available - this example is 
selected for demonstration purposes to allow a simple and clearly defined convergence criterion). 

For each of the parameter sets, 1,000 calculation sets were performed using RSS, HLHS, and 
SRS. Each calculation set begins with 20 samples and is extended until the convergence criterion 
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Excess Kurtosis 


(a) 


(b) 


Figure 11: (a) Number of sample required to obtain 95% convergence using SRS, HLHS, and RSS as a function 
of response kurtosis. (b) Reduction in sample size for HLHS and RSS compared to SRS as a function of response 
kurtosis. 


is satisfied. The number of samples required for convergence is recorded for each set. The results 


of this study are summarized in Figure 11 Figure 11(a) shows the number of samples required to 
achieve convergence for 95% of the calculation sets using SRS, HLHS, and RSS as a function of 
response excess kurtosis. In general, a larger number of samples are required for convergence as 
the kurtosis increases although the proposed RSS method consistently outperforms both SRS and 
HLHS in this regard. Figure [l 1 (b)| shows the relative reduction in sample size over SRS for both 
HLHS and RSS defined by: 

Ng — IV h! R 

h/r ( 2 7 ) 


N S 


where Ng, N H / R are the number of samples required for 95% convergence (plotted on the left) 
for SRS and HLHS/RSS respectively. By this measure, the RSS method affords upwards of a 
90% reduction in sample size when compared to SRS for responses with low kurtosis. Initially, the 
effectiveness drops precipitously as the kurtosis increases. After kurtosis ~ 2, it stabilizes while still 
providing significant sample size reduction. Notice also that RSS consistently produces a 10-20% 
greater reduction in sample size than HLHS. A complete tabular summary of this study is provided 
in the Appendix. 

Next, we consider the convergence of the empirical CDFs. The “true” CDFs Dy(y) for each 
case A-J are determined by Monte Carlo simulation with 100,000 stratified samples. Convergence 
to Dy(d ) is measured for the first 10,000 samples using the area validation metric (a.k.a. the 
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Minkowski L\ norm) as |46j: 


/ °0 

\D Y {y)-b^\y)\dy. (28) 

-OO 

where Dy (y) denotes the empirical CDF computed from l samples. Sample plots showing the 
convergence of 5(1 ) using SRS, HLHS, and RSS are provided in Figure [12} While HLHS and RSS 
provide comparable rates of convergence, HLHS allows only those sample sizes shown with an ‘x’; 
highlighting the need to add many samples during an HLHS sample size extension. RSS meanwhile 
affords extension by as many or as few samples as desirable for the application. 



Number of Samples 



10 1 10 2 10 3 10 4 


Number of Samples 


Distribution E 




Figure 12: Convergence of the empirical CDF to the “true” CDF for select cases using SRS, HLHS, and RSS. 


6.2. Structural response to underwater shock 


Growing emphasis is being placed on uncertainty quantification in computational physical mod¬ 
eling and simulation where it is desirable to evaluate the response variability for large and complex 
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physical systems. We consider the case of the multi-physical response of a floating structure to un¬ 
derwater shock. These calculations involve large Lagrangian structural finite element models (with 
tens of thousands of elements) coupled with extremely large Eulerian fluid models (with tens to 
hundreds of millions of fluid elements) as depicted in Figure [Tifjyi). Such calculations often require 
several days to weeks of computation time on massively parallel computers. An accompanying 
paper on the validation of computational models for this problem was recently published by the 
authors m • 


m 



y Lagrangian 
/ Structural Model 


Eulerian Fluid Mesh 


(a.) 


* 



Figure 13: (a.) Schematic of coupled fluid structure model for underwater shock response of a surface vessel, (b.) 
Simplified model used for demonstration purposes. 


To demonstrate the tractability of uncertainty quantification for these problems, it was necessary 
to show that statistical convergence could be accomplished for the response in a small number of 
samples. To do this, a simplified representation of the system was constructed as shown in Figure 
[l3|(b.). The simplified model is an uncoupled one dimensional, two degree of freedom damped 
system consisting of a mass M\ representing a primary structure supporting a much larger piece 
of equipment with M 2 = 5M 1 . Simplified representations of this form have been used for decades 
to study fluid-structure interaction in near-surface underwater shock. Early work by Bleich and 
Sandler [I5J presented the effects of bulk surface cavitation on the velocity response of a floating 
structure subject to underwater shock in a bilinear cavitating fluid. This was followed by the 
development of a numerically uncoupled treatment of the surface shock response problem in ID 
that accounts for the effects of cavitation by DiMaggio et al. |29]. I n a few more recent papers 
describing a numerical technique for capturing the effects of cavitation in finite elements models, 
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Sprague and Geers EDI EH, with corrections by Stultz and Daddazio [52], presented the two mass 
problem (neglecting the effects of structural damping) as it is studied here. 

Unlike the previous authors, we are not specifically concerned with the effects of bulk cavitation 
and the complex physics that drive the response following the initial kickoff. Instead, we are 
interested in capturing the variability in peak velocity of the mass M 2 given uncertainties in the 
model using a minimal number of samples. This is motivated by the need to verify that simulation- 
based uncertainty quantification is a tractable option for the more complex multi-physics model. By 
considering the sensitivity of the model to various parameters and assessing realistic uncertainties 
associated with this problem, we identify two uncertain quantities. First, damping of real structures 
is notoriously difficult to ascertain and model with any accuracy. This is further complicated by 
the dissipation of energy resulting from fluid-structure interaction. Linear viscous damping (£^) 
is assumed and considered to possess a truncated (non-negative) normal distribution with mean 
/id = 0.025 and standard deviation = 0.01. Second, noticeable variabilities are known to arise 
in the peak pressure produced from underwater detonation of certain explosives. Here, we consider 
charge weight to be a normal random variable with mean /i c = 117 lb. TNT and standard deviation 
a c = 1.17 lb. TNT to account for this variability in peak pressure. 

6.2.1. Modified bootstrap convergence evaluation 

To evaluate statistical convergence, a bootstrap procedure is used that enables us to estimate 
approximate confidence bounds for response statistics of interest. The traditional bootstrap pro¬ 
cedure, developed by Efron [53], performs a random resampling (with replacement) of a dataset 
to obtain a large number of surrogate sample sets. Statistical analysis of these so-called bootstrap 
samples provides an estimate of the variability of the computed statistic. The traditional boot¬ 
strap method resamples from the original dataset with each sample possessing equal probability 
of occurrence. That is, the probability of selecting sample Si from a set of N samples is given 
by P[Si] = 1/N. However, in stratified samples, and more specifically the RSS procedure devel¬ 
oped herein, the samples possess associated weights that are not necessarily equal. To account for 
this unequal weighting (imbalanced stratified design), a modified bootstrap method is proposed as 
follows: 

1. Find the least common denominator D for all sample weights and create a modified dataset 
of size D by duplicating samples such that each sample in the modified set has equal weight. 
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2. Resample a set of N points (bootstrap sample) with replacement from the modified set. 

3. Compute the statistic of interest (bootstrap statistic) from the bootstrap sample. 

4. Repeat a large number of times to obtain an empirical distribution for the bootstrap statistic. 


The proposed resampling process accounts for the sample weights but, it should be emphasized 
that the bootstrap datasets produced from this method are not stratified samples. In any given 
bootstrap resample, it is likely that certain strata will not be represented. The consequence of this 
loss of stratification is a potential increase in the variance of the bootstrap statistics. 

This new bootstrap method provides an admittedly imperfect monitor of convergence. The 
shortcomings of the proposed monitor motivate the need for further research to establish a more 
effective and general monitor. 


6.2.2. Results 

The bootstrap procedure outlined in the previous section is used to determine confidence inter¬ 
vals for the response statistics from the empirical CDF of the bootstrap statistics. We are interested 
in evaluating statistics of the peak velocity of mass M 2 denoted by V 2 . Convergence is defined such 
that the bootstrap 95% confidence bounds for each statistic of interest are sufficiently narrow. 
Denoting a computed response statistic by T with bootstrap 95% confidence intervals given by 
Tg 5 and Tg 5 corresponding to the lower and upper bound respectively, the convergence criterion is 
expressed as follows: 


rpu rpl 

^ 95 1 95 


T 


< e t h 


(29) 


Note that the convergence criterion based on confidence intervals from an empirical CDF was 
selected in order to account for skewness of the likely non-Gaussian bootstrap CDF. A convergence 
criterion based on variance of the bootstrap statistics, for example, would not account for this 
reality. 

Convergence was compared for RSS and SRS. Considering the performance comparison present 
previously and the need for maximal flexibility in minimizing sample size, HLHS was not considered 
a good candidate for these analyses. Initially 90,000 analyses were performed using SRS and the 
convergence evaluated. To achieve superior convergence in all statistics, far fewer calculations were 
required using RSS. In all, 4,000 RSS analyses were conducted although less than this were necessary 
to match the convergence or SRS. Plots of the first four moments (with bootstrap confidence 


29 



intervals) as a function of sample number for peak velocity V 2 are provided in Figure [Id] using both 
RSS (left) and SRS (right). 
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Figure 14: Evolution of response statistics with added samples for peak velocity of mass M 2 using SRS (left) and 
RSS (right) showing 95% bootstrap confidence intervals for 4,000 samples. 
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Table 4: Number of samples required to achieve various levels of convergence, -’s indicate that this level of convergence 
was not obtained in the number of samples performed (4,000 for RSS and 90,000 for SRS). 


Sy 2 

15% 

5% 

2.5% 

1% 

0.5% 

0.1% 


SRS 

< 20 

40 

140 

840 

3,480 

84,900 

RSS 

< 20 

< 20 

< 20 

60 

140 

300 

a 

SRS 

460 

4,500 

18,800 

- 

- 

- 

RSS 

120 

280 

460 

1,760 

3,760 

- 

7i 

SRS 

10,750 

73,350 

- 

- 

- 

- 

RSS 

1,760 

2,240 

- 

- 

- 

- 

72 

SRS 

69,700 

- 

- 

- 

- 

- 

RSS 

2,920 

- 

- 

- 

- 

- 


By establishing convergence thresholds, we determine the number of samples required to produce 
sample statistics that are considered adequate. Table [4] shows the number of samples required for 
various convergence criteria of the first four statistical moments for both RSS and SRS. It is 
noted that the mean value converges extremely rapidly with a cutoff threshold of 0.1% requiring 
only 300 samples for RSS compared to nearly 85,000 for SRS. Furthermore, a reasonably small 
error can be gained in standard deviation with only hundreds of samples. Skewness and Kurtosis 
meanwhile require thousands of samples to converge even using RSS. Nonetheless, the improvement 
in convergence over SRS is an order of magnitude or better. 

These data are encouraging when considering the prospects for an uncertainty quantification 
study for the large-scale multi-physics problem, particularly if the study is interested in evaluating 
low-order statistics (mean and standard deviation) of various response quantities. As indicated 


by the plots in Figure 14, the method converges very rapidly for these statistics even when the 
response quantity of interest (here peak velocity of mass M 2 ) deviates appreciably from Gaussian. 
As a result, this large-scale multi-physics study was performed but details of this study are beyond 
the scope of this paper. 


7. Conclusions 

In this work, we promote an adaptive approach to Monte Carlo-based UQ rooted in stratified 
sampling. To motivate the advantages of SS for certain classes of problems the space-filling, orthog- 
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onality, and projective properties of SS are studied and compared with simple random sampling 
(SRS) and Latin hypercube sampling (LHS). It is shown that SS provides superior properties to 
many low-to-moderate dimensional problems - especially when strong variable interactions occur. 
To enable the adaptive approach, a new sample size extension methodology - Refined Stratified 
Sampling (RSS) - is proposed that adds samples sequentially by dividing the existing strata. RSS 
is proven to reduce variance compared to existing sample size extension methods for SS that add 
samples to existing strata and is shown to possess comparable or improved convergence when com¬ 
pared to existing extension methods for LHS while affording maximal flexibility. Using RSS, sample 
can be added one-at-a-time while sample size grows exponentially using hierarchical LHS. Opti¬ 
mality of the stratum refinement is examined and its extension to high dimensional applications 
discussed. Several future research opportunities along these lines are identified that may facilitate 
significant breakthroughs in sample-based UQ if achieved. In particular, the potential ability to 
adaptively identify a partitioning (stratification) of the input space that optimizes the partition¬ 
ing of the output space could serve to minimize sample size dramatically compared to the current 
state-of-the-art. Additionally, the RSS method, when combined with new methods for decomposing 
the sample space into low-dimensional subspaces will allow its extension to very high dimensional 
problems. 

Two examples are presented. The first involved the forward propagation of uncertainty through 
a low-dimensional stochastic function where convergence can be explicitly evaluated at each sample 
size extension. The second involved a real multi-physics computational model possessing uncertain 
parameters for the response of a floating structure subject to an underwater shock. In the sec¬ 
ond example, a new bootstrap procedure was proposed for resampling from a stratified design to 
compute confidence intervals for response statistics. 


Appendix 


Example 1 analyzes the number of simulations required to achieve convergence for the 3- 


dimensional stochastic function given in Eq. (25) with 10 different sets of input distributions 


using SRS, LHS, and RSS. Table [5] provides the results of this study for all three sample size exten¬ 
sion methods and each set of input distributions defined in Table [3| The columns show the number 
of samples required to achieve different proportions of converged calculations. For example, 10% 
of the 1,000 (or 100) calculation sets converged within 72, 40, and 42 samples for Distribution A 
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using SRS, HLHS, and RSS respectively. In general, the RSS method shows a drastic reduction 
(often greater than an order of magnitude) in the number of samples needed for convergence when 
compared to SRS and even affords a significant savings when compared to HLHS. In particular, 
it is apparent that oversampling is necessary with HLHS given its exponentially increasing sample 
size. 
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